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We have applied the Melnikov criterion to examine a global homoclinic bifurcation and transition to chaos in a case of a 
double well dynamical system with a nonlinear fractional damping term and external excitation. The usual double well 
Duffing potential having a negative square term and positive quartic term has been generalized to a double well potential 
with a negative square term and a positive one with an arbitrary real exponent q > 2. We have also used a fractional 
damping term with an arbitrary power p applied to velocity which enables one to cover a wide range of realistic damping 
t factors: from dry friction p — *• to turbulent resistance phenomena p = 2. Using perturbation methods we have found a 

critical forcing amplitude fi c above which the system may behave chaotically. Our results show that the vibrating system 
is less stable in transition to chaos for smaller p satisfying an exponential scaling low. The critical amplitude fi c as an 
f"*^ , exponential function of p. The analytical results have been illustrated by numerical simulations using standard nonlinear 

tools such as Poincare maps and the maximal Lyapunov exponent. As usual for chosen system parameters we have 
identified a chaotic motion above the critical Melnikov amplitude /i c 

Ph ' 
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> \ 1 Introduction 

A nonlinear oscillator with single or double well potentials of the Duffing type and linear damping is one of the simplest 



systems leading to chaotic motion studied by (T]|2||3j[4|. The problem of its nonlinear vibrations has attracted researchers 
from various fields of research across natural science and physics 1 4 5 6 1, mathematics 1 8 1 mechanical engineering 1 10 11 
, I14II12|[T31 : and finally electrical engineering (l]|2||3. This system, for a negative linear part of stiffness, shows homoclinic 
orbits, and the transition to chaotic vibration can be treated analytically via the Melnikov method |7|. Such a treatment has 
been already performed successfully to selected problems with various potentials |8 9 14|. Vibrations of a single Duffing 
oscillator have got a large bibliography fl]l2|3]|l|5]©|l|l[IO|[ni[I2[IlE last decade coupled 

Duffing oscillators 1 19, 20, 21 22 23 1 with numerous modifications to potential and forcing parts have been studied. On 
the other hand the problem of nonlinear damping in chaotically vibrating system has not been discussed in detail. Some 
insight into this problem can also be found in the context of self excitation effects 

1 1 51 1 1 91 1221 1231 1251 1241 1261 . and dry friction effects 1 331 1341 1351 1361 In the paper by Trueba et al. 1161 . the systematic 
discussion on square and cubic damping effects on global homoclinic bifurcations in the Duffing system has been given. 
Recently Trueba et al. 1171 and Borowiec et al. 1 18 1 have analyzed a single degree of freedom nonlinear oscillator with 
the Duffing potential and fractional damping. Different aspects fractionally damped systems have been studied recently 
by Mickens et al, Gottlieb, and Mickens |27, 28, 29 1. On the other hand Maia et al. and Padovan and Sawicki 1 30 
1321 1311 analyzed similar systems where fractional damping have been introduced in different way through a fractional 
derivative. Awrejcewicz and Holicke |37| and more recently Awrejcewicz and Pyryev |38| applied Melnikov's method 
in the presence of dry friction for a stick-slip oscillator. More general introduction to the problem of non-smooth or 
discontinuous mechanical systems can be found in 1391 1401 . 

In the present paper we revisit this problem looking for a global homoclinic bifurcation and transition to chaotic vi- 
brations in a system described by a more general double well potential where its usual positive quartic term has been 
generalized to term with an arbitrary real exponent grater than 2. Below we would also apply a nonlinear damping term 
with a fractional exponent covering the gap between viscous, dry friction and turbulent damping phenomena. 
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The equation of motion has the following form: 

x + ax \x\ p 1 + Sx + 7Sgn(a;)|a;| 9 ~ 1 = ficosuit, 
where x is displacement and x velocity, respectively, while the external force F x : 



(1) 



F x = -Sx - jsgn(x)\x\ q \ 
and corresponding potential V(x) (Fig. ^) is defined as: 

2 q 



(2) 



(3) 



where q > 2 is a real number. In spite of the definition V(x) (Eq. [3} in terms of absolute value |x| it is still a function of 
C 2 class if only q > 2 (see Appendix A). 

The non-linear damping term is defined by the exponent p: 



dpt(i) = ax \x\ p 



(4) 



In Fig. \\]p we have plotted the above function versus velocity (v — x) for few values of p. Note that, the case p — + (see 
p = 0.1 in Fig. [\]p for a relatively small velocity) mimics the dry friction phenomenon [33 34) . 
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Fig. 1 External potential V(x) = -| — h 71 1 (Eq. 3) for 5 = —2 for a few values of g (g = 7 > 2 in Fig. la), Damping term for 
various p (Fig. lb). 



2 Melnikov Analysis 

We start our analysis with the unperturbed Hamiltonian: H° 

H" = ^ + V(x). (5) 

Note that for our choice of potentials 5 = —2 and 7 = q (Fig. [^) V(x) has the three nodal points (x = —1, 0, 1) where 
the middle one (x = 0) corresponds to the local peak at the saddle point. The existence of this point with a horizontal 
tangent enables occurrence of homoclinic bifurcations. This includes transitions from regular to chaotic solutions. To study 
the effects of damping and excitation on the saddle point bifurcations, we apply small perturbations around the homoclinic 
orbits. Our strategy is to use a small parameter e to the Eq. Qwith perturbation terms. Uncoupling Eq. [Qinto two differential 
equations of the first order we obtain 

x = v (6) 
v = — eav\v\ p 1 — 5x — 7Sgn(a;)|a;|' J_1 + eft, cos uit, 

where ea = a and eft — /1, respectively. 
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Fig. 2 Left hand side homoclinic orbits for unperturbed Hamiltonian (Eq. Note, in our case the potential has reflection symmetry 
over 0-y axis so the orbits appear in pairs for corresponding regions x > and x < 0. 



At the saddle point x — 0, for an unperturbed system (Fig. [^), the system velocity reaches zero v — (for infinite time 
t = ±oo) so the total energy has only its potential part which has been gauged out to zero too. Thus transforming Eqs. [3] 
and[3]for a nodal energy (E = 0) and for 5 < 0, 7 > we get the following expression for velocity: 



Sx 2 , y\x\ q 
"~2 V 



Performing integration over x we get 

d.r 



t-tn 



± 



-<y- 



2-j\x\i- 2 

1 



(7) 



(8) 



where to represents here a time like integration constant. 

Integration in Eq.|8]has been performed analytically. For q > 2, one can write x* as: 



x*(t-t ) = ± 



-5q \ " 
^7 



1 



coshi- 2 



^V=7(< - to) 



(9) 



The corresponding velocity v* reads: 

v* = v* (t — to) = ^fVS 



-Sq\ 
~2T 



tanh 






cosh?- 2 


^>/=7(* - to) 



(10) 



Due to the reflection symmetry of potential V(—x) = V(x) (Eq. there are two symmetric solutions for unperturbed 
homoclinic orbits with '+' and '-' signs in Eqs. 191 101 A family of right hand side homoclinic orbits (x* , v* ) has been plotted 
in Fig. I2 In unperturbed case both stable and unstable manifolds (Poincare sections of the orbits are usually denoted by 
W s and W u ) can be identified with the orbits discussed above while perturbations would influence them in a different way 
1121 . Existence of cross-sections of between W s and W u manifolds signals Smale's horseshoe scenario of transition to 
chaos. 

The distance d (Fig. between them can be estimated by the Melnikov function M (to): 



M(t ) 



+ OO 



h(x*,v*) f\g{x*,v*)dt 

where the corresponding differential form h means the gradient of unperturbed Hamiltonian (Eq. 

h = (8x* + 7Sgn(a;)|a;|' 3 " 1 ) dx + v*dv, 



(11) 



(12) 
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Fig. 3 Schematic plot of stable and unstable manifolds (W s and W u ) of perturbed system Eq. [6] d denotes the distance between 
manifolds given by Melnikov function M (to) Eg. II II 



while g is a perturbation form (Eq. |5} to the same Hamiltonian: 

g = (^flcosut — av* \v*\ p ^ dx. 

All differential forms are defined on homoclinic orbits (x, v) = (x* , v*) (Eas. 191101 
Thus the Melnikov function M (to): 

/+oo 
v*(t) [p,cos(u(t + t ))-av*{t) |u*(i)| p_1 J dt 



(13) 



+ 00 



■sina;io / v*(t)fl sin wtdt 

J — oo 

sin(a;to)Mi — a/ 2 , 



+oo 



av^^lv^t^dt 



(14) 



where I\ and I2 are integrals to be evaluated. sinw<o appears because of the odd parity of the function v*(t) under the 
above integral where 



cos(u(t + to)) = cos(ujt) cos(wto) — sin(ojt) sin(wio)- 



(15) 



Thus a condition for a global homoclinic transition, corresponding to a hors-shoe type of stable and unstable manifolds 
cross-section (Fig. 0, can be written as: 



\jM(t o )=0 and 

to 



dM(t Q ) 
dtn 



^0. 



(16) 



For a perturbed system the above constraint together with the explicit form of Melnikov function Eq. [Hlgives the critical 
amplitude /i c : 



tW_ 

a 



where 1\ and I2 are corresponding integrals given in Eq.^J In case of 1\ we have the following integral 



h 



-Sq \ " 
~2j 



-00 tanh (^^y^f 



cosh-j- 



sin(wi) dt 



(17) 



(18) 



to be evaluated numerically in general but for some cases can be easily performed numerically (see Appendixes B and C) 
while, in analogy to 116111711181 . 12 can be expressed as 



h = (-6) 



I ,+oc sinh p+1 [i^l^—5t 
2 7 J loo cosh^ 1 (^V^ 



p+: 

■ 1 / — 8q \ q -' 



dt 



27 



'P+D 



B 



p+2 p+1 



2 '(q-2)J> 



(19) 
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Fig. 4 Critical amplitude [ic/at versus frequency u> for few values of p (p = 0.1 0.5, 1.0, 2.0) and different q (q = 2.5 in Fig. [4jt, 
5 = 3.0 in Fig. EJ>, q = 3.5 in Fig. gj;, p = 4.0 in Fig. |4jl. ln(/i c /a) versus the exponent p for three values of to and q = 2.5 and 3.5 
(Fig. g^). Dependence of the slope k (k = tan</> in Fig. g£) on the exponent q - Fig. |4f (T corresponds to present calculations for 
ui = 1 while '2' is a fitting trail k = l/(g/1.3 + 1)). 



where B(r, s) is the Euler Beta function dependent of arbitrary complex arguments with real parts (Rc r > and Rc s > 0) 
defined as 

r(r + s) 

while r(r) denotes the Euler Gamma function: 



/•oo 

r(z) = / , 

Jo 



e~ s s z - L ds for Re z > 0. (21) 

In Fig. |4^-d we plotted the results of Melnikov analysis for a critical amplitude fi c /a for few values of p (p = 0.1 0.5, 
1.0, 2.0) and different q (q = 2.5 in Fig. gt, q = 3.0 in Fig. gj), g = 3.5 in Fig. gk p = 4.0 in Fig. gjl) For /it > ^ c the 
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Fig. 5 Maximal Lyapunov exponent Ai versus fi (Fig.|5p), bifurcation diagram together with size of attractor x max and x m in versus fi 
(Fig.|5j3)forw = 1.1. 



system can transit to chaotic vibrations. Note, in spite of some quantitative changes all four figures (Fig. |4^-d) have similar 
shape and the sequence of corresponding curves with particular exponents p = 0.1 0.5, 1.0, 2.0 is preserved for any u and 
q. This gives us a conviction that p may play some independent role. In fact plotting ln(/i c /a) in Fig. versus p and we 
have got straight lines with characteristic slope independent on uj but changing with q. Dependence of the slope n: 



tan < 



(22) 



defined in Fig. versus q has been plotted in Fig. |4f for u> = 1. Note, the curve '1' corresponds to present calculations 
for while '2' is a fitting curve: 



1 



1.3 



(23) 



The above scaling is not a surprise taking into account the structure of M (to) (Ea. ll2l >. In this expression the exponent p is 
entering to the second integral independent of u>. On can also look into the analytic formulae for n c in cases of q = 4 and 
3 in the Appendix B (Eq. IB.5IB.8t where the p appears as an exponent. 



3 Results of Numerical Simulations 

To illustrate the dynamical behaviour of the system it is necessary to simulate the proper equations. Here we have used 
the Runge-Kutta method of the forth order and Wolf algorithm |41 1 to identify the chaotic motion. In our numerical code 
we started calculations from the same initial conditions (xo,t>o) = (0.45,0.1) for any new examined value of fj,. The 
system parameters S — —2 and 7 = q have be chosen the same as for analytic calculations. We have performed numerical 
calculations for different choices of system parameters: a, ui, p and q, hut here, for or technical reasons, we limited our 
discussion to a = 0.1, to = 1.1, p = 0.5 and q = 2.5. 

In Fig. [5^ we have plotted the maximal Lyapunov exponent versus external forcing amplitude /i. Here one can clearly 
see points of Ai sing changes. For /1 G [0.23, 0.27] and [0.33, 0.38] we have got Ai > indicating chaotic vibrations. In 
Fig. |5j3 we have plotted the corresponding bifurcation diagram. Thus a black region means chaotic motion. This result, as 
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Fig. 6 Critical amplitude fx c /a versus uj the dashed line corresponds to lj — 1.1 (Fig. |6^). Phase portrait and Poincare maps for 
lu — 1.1, q = 0.1 and two different [i (p = 0.05 in Fig.|6j) while fi — 0.24 in Fig. [(J). The results have been obtained forp = 0.5 and 
q = 2.5. 



well as others, calculated for different sets of system parameters, is consistent with the Melnikov results. For comparison 
we have plotted the Melnikov curve again (Fig. |6^ with two trial points fx = 0.05 and fx = 0.24 (for a — 0.1). There is no 
doubt that Fig. |6ji shows the regular synchronized motion represented by a single loop on a phase portrait and a singular 
point on Poincare stroboscopic map. On the other hand Fig. |6j; shows clearly a strange attractor of chaotic vibration with 
complex structure of the Poincare map. 



4 Summary and Conclusions 

We have examined criteria for transition to chaotic vibrations in the double well system with a damping term dpt(w) = 
v\v | p_1 described by a fractional exponent p and nonlinear potential with negative square term (related negative stiffness) 
and a positive term with higher exponent \x\ q where q > 2. In spite of non-smoothness of corresponding vector-fields 
(h and g - Eqs. ^2 an dEl respectively) it has been proven in the Appendix A that extra terms to the Melnikov integral 
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1 40 1 projected out. Thus the critical value of excitation amplitude fj, above which the system vibrate chaotically has been 
estimated, in by means of the Melnikov theory |7|. For some selected values of the exponent q (q —4,3, 8/3, 2/5) it was 
possible to derive a final formula for \i c while for other cases one of Mielnikov' integrals has be calculated numerically. 

The analytical results have been confirmed by simulations. In this approach we used standard methods of analysis as 
Poincare maps, bifurcation diagrams and Lyapunov exponent. 

The Melnikov method, is sensitive to a global homoclinic bifurcation and gives a necessary condition for excitation 
amplitude fi = \i c \ system in its transition to chaos 1 8 9 1. On the other hand the largest Lyapunov exponent |41 1, measuring 
the local exponential divergences of particular phase portrait trajectories gives a sufficient condition /.i = /i C 2 for this 
transition which has obviously a higher value of the excitation amplitude fi = fx C 2 > Mci ■ 

Above the Melnikov transition predictions (/i > /i c i) we have obtained transient chaotic vibrations !9l ll0lfTTlll2lll8l 
as we expected drifting to a regular steady state away the fractal attraction regions separation boundary. This is typical 
behaviour of the system which undergo global homoclinic bifurcation. 
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Appendix A 

Starting with the perturbation equation (Eq. |6jl we write it in a two element vector form 

q = h + eg, (A.l) 

where 

q = [x, v] 

h = [v,-Sx- 1 sgn(x)\x\ q - 1 ] (A.2) 
g = [0, — av H P_1 + flcoscot]. 
On the other hand the homoclinic orbit 

q*(t-t ) = [x*(t-t ),V*(t-t )], (A3) 

where to is usually defined by simple zero of Melnikov integral |7|. In the limit of extreme time t — ► ±oo the system state 
[x, v] reaches a saddle point [x, v] = [0, 0] (see Figs. [TJi, [5J|. Consequently, in the aim to examine the Melnikov criterion 
for chaos appearance, the vector fields h and g are defined on the homoclinic orbit (Fig. |2j as: 

h(q*) = [v*(t-t ),-Sx*{t-t )-'yBea{x*{t-h)\x*{t-t )\ q ~ 1 ] (A.4) 

and 

g(q*,t) = [Q,-av*(t-t ) \v* {t - to)]^ 1 + fl cos ut], (A.5) 
The perturbed stable and unstable manifolds W s and W u read |40| 

q"< s (i, t ) = q*(t - to) + eq?" (Mo) + 0(e 2 ) (A.6) 

respectively. 

Note the perturbation correction to the homoclinic orbit qi' s (t, to) in the above expression (Eq. IA.6> should be found 
by solving the following linear differential equation about the examined time t (or a system state q = q* (t — to)): 



q?*'(*.*o) = ( ir - ir ) q r(*.*o) + g(q*(*-*o),*) (A.7> 

av ax ) | q=q , (t _ io) 



Note that the above vector fields h(x, v) and g(a;, v, t) are not C 2 functions. Namely h is of C 2 only if q > 2 and of C 1 
if 1 < q < 2. Similarly g is of C 2 only if p = 1 or p > 2 and of C if < p < 1 and C 1 1 < p < 2, respectively. In case 
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of 1 < q < 2 the line x = separates the whole phase space (x, v) into two parts where h(a;, v) is enough smooth (C 2 ). 
The same can be applied to the line v = and g(x, v, t) as a possible set for < p < 1 and 1 < p < 2. This line crosses 
manifold at [xd, vj\ for the specific time t = td such that Xd — x(td) and Vd — v{td) (Fig. IA.ll . 

According to Kunze and Kiipper [40 1 the (C 2 ) space separation includes additional terms to the Melnikov integral. 

Thus the Mielnikov function M(io): 



M(t )=Mo(to) + ^(h ± - + (q*(^))-q^ + (io+^ ! to)-h ± -(q*(^)).q^-(t + ^,i ) (A.8) 

+ h+"(q*(t+)) • qj'-(i + t+,t ) - h+"(q*(t+)) ■ q* ,+ (i + t+,t )) , 

where q*' , q"' are stable and unstable manifold perturbation solutions (Ea. lA.7> for t in the vicinity of td but t > td for 
'+' sign and t < td for ' — ' sign, respectively. 
Mo (to) is defined as for smooth vector fields: 

/oo 
h X (q*(t-t )-g(q*(i + io),i)crf (A.9) 
-oo 

and rr 1 = [— hz, hi]. 

Once first discontinuity are identified for x = 0, td — ► ±oo, one have to examine <f>f(t,to) — 5i' s ' (t, to) and 
<f>t(t,t ) = q? s < ± (t,t ): 



^ = (l + ^ + 7(«-l)^- 2 )^i for x>0 (A1Q) 

h = (1 + 5 + j{q- l)xi- 2 )<p 2 -avlvl?- 1 +ficosut 

and 

^ 1 = (l + 5 + 7 ( g -l)(-^- 2 )0 1 
/. 2 = (l + <5 + 7( g -l)(-a;)«- 2 )0 2 -(5t>| U |P- 1 +/icoswt ' k ; 

Note substituting x = for q > 2 to Eqs. IA.lOl and lA.1 ll we get the same equations for <f>i/ 2 and consequently the same 
expressions. This means automatically no extra terms to the Melnikov integral (Eq. IA.8t caused by the x — discontinuity. 
Interestingly q < 2 would lead to a different result but for such case the is no homoclinic orbit in the unperturbed system 
described by H° (Eqs. EE}. 

Let us non focus on v = discontinuity. In this case 

j>i = (1+^+7(5- i)k*r 2 )^ for v>0 (Am 

<t> 2 = {l + 5 + 1 (q-l)\x*\i~ 2 )<j }2 -a(v*y + pcosujt V ' ' 



fc = (1 + ^ + 7(5- l)|a1*- 2 )0i for V<Q (Am 

= {l + 5 + 1 {q-l)\x*\*- 2 )<j )2 +a{v*y + ficosujt ^ K ' 
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Note, excluding natural odd numbers for the p exponent, the above equations (Eas. IA. 12l and lA~T3t are usually different 
for any other p > 0. However both solutions [4>i , </>2~] and (f)^] have to be projected into 

K-=o = [-ft2,/ii]|„=o = [Sx*(t d - t Q ) + 7 sgn(x*(t - t )\x*{t d - t )| 9_1 ,0] (A.14) 

and the differences in solutions in 02 an d 4>2 me effectively projected out. Interestingly this is also valid for p = (a dry 
friction case). 

Finally for q > 2 and p > the Mielnikov function M(irj) can be treated as a 

M(to)=Mo(to). (A.15) 



Appendix B 



In this appendix we show how to get homoclinic orbits and analytically for some specific cases of exponent q: q = 4, 3, 

2.67 and 2.5. 

In case of q = 4 we follow works by Trueba et al. (T6 1 and Borowiec et al. Jl 8 1 (and Eqs. 191101 



x* = x*(t — to) = ±4/ 



1 



v* = v*(t-t ) 



7 cosh (-\/— S(t - t )) 
-u fix tanh (V^O - to)) 



(B.l) 



7 cosh - h)) 

where '+' and '— ' signs are related to left- and right-hand side orbits, respectively, to is a time like integration constant. 
On the other hand for q = 3 we have 

35 



x* = x*(t — to) = T 



v* =v*(t-t ) = T 



1 



2 7 cosh 2 h 



tanh 



S(t-tg) 



2 7 cosh 2 (- 



(B.2) 



Consequently for q = 2| = 8/3 ~ 2.67 we have 



45 

X*=X*(t-to) = ±(-3- 



.3/2 



cosh 



i;* =«*(*- to) 

And for q = 2.5 

2;* = x*(t — to) 



-46 
"37 



3/2 



tanh 



x -<5 



/=5(t-tp) \ 

3 J 



cosh' 



(B.3) 



v* =«*(*- to) = ± — V-<5 



4 7 y C osh 4 ( ^- to) ) 

( 



45 



3 



tanh 1 



-g(t-to) 



cosh' 



^ y/=?(«-to) ^ 



(B.4) 



The results for a Melnikov integral can be easily found in the above cases. Evaluating the corresponding integral (Eq. 
II It after some algebra the last condition (Eq. I16> yields to a critical value of excitation amplitude /i c . Thus for q = 4 
HEiriliniD we have: 



2 P/2(_,5)P+l/2 /p+2p+l\ / 

Mc = « — o I — - — , — - — I cosh 



p/2 



2V = 5 



(B.5) 
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while in case of a = 3 l42l ?.l43l: 

3p(-,5) 3 p/ 2 +2 fp + 2 \ . , 



2 f(p+i)(_ 5) n P /io-2/5 /pinip+m / 3ttw 

= a , , — — : ; - _ — d I — - — , I cosh 



forg = 8/3: 



3*(p +1 )tt(9cj 2 - <5) W7 3p/5+9/io V 2 2 
and finally for q = 5/2: 



5 2 P( _ (5) 5p/2 + 3/2 fp+2 A /2tTW 

Mc = a ^,.,, . 9 K .,.,„„ B ( — ^— , 2(p + 1) smh 



"2 4 p+ 3 tt(4w 2 -5)cj 2 7 2p V 2 



-J 



2V^y ' 



(B.6) 



(B.7) 



(B.8) 



Appendix C 



The integral 1\ can be evaluated analytically in some specific cases of exponents q corresponding to homoclinic orbits 
Eqs. IB . 1 IB .41 numbered by the corresponding power index to applied to hyperbolic cos function in the denominators. Let 
us consider integrals ij for given to = 1,2,3 and 4 related to various q exponents q = 4, 3, 8/3 and 5/2, respectively. To 
better clarity we will use new notation I\ — » 1\ (m) for given m: 



+ °° n- ■ , „ /" +0 ° tanh(r) . , 

v (tj/zsinurfdi = C m / m . . sin(a> m Tjdr 



x cosh m (r) 
f 00 cos(u m -r) w m C„ 



h(m) 



to cosh m (r) 
while constants C m and u) ra tire defined as follows: 

m/2 



(CI) 



dr 



TO 



Cn 



-s 



-(to + 1)<5 



(C2) 



TO7 y v— <5 

Evaluating the integral J m (uJ m ) iC.H . for positive integer m, twice by parts we have got the following recurrence identity 



Jm+2{u r n+2) 



to(to + 1) 



for 



to = 1, 2, 3, . 



(C3) 



Thus only J\ and J-2 need to be calculated. Below we evaluate them on the complex plane by summing corresponding 
residue. 



y 


complex plane 
z=x+iy 


y 


/ k=3 i 






/ 2 






/ 1 , 






1 1 




V 


-1 1 


X 




-2 < 







Fig. C.l Deformed contour integration schema and imaginary poles. 



J k=i 



(C.4) 
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where 



Res [/(«), z fc ] 



lim 



(to — 1)! z^zk dz 



T [(z-z fe r/(z)]. 



(C5) 



for to — 1 or 2, in our case. 

The examined function f(z) is defined as: 



/(*) 



■ exp (vjj m z). 



(exp(z) + exp(-z)) r 
Note that on the real axis (Fig. IC.H Rc z = t it can be written as 

cos(cj m r) 

Im/(T) = rfv • 

The multiplicity of each pole of the complex function f(z) (Eq. IC.6I is given by 



z k = (Z + i 



for 



fc = 1,2,3,. 



(C6) 



(C7) 



(C8) 



Note J m (Eq. IC.U can be easily determined for to = 1 or 2. Namely, after summation of all poles in the upper half-plane 
(Fig. IC.H . we get for to = 1 



while for to = 2 we obtain 



Ji = 



+ 00 



coshr cosh(2|i) 



7TCJ 2 



cos(o;2t) 
cosh 2 r ~ sinh(^) 



(C9) 



(CIO) 



On the other hand, in case of to = 3 and to = 4 (and also for any larger to), we can use the recurrence relation (Eq. 
ED: 



Consequently using Eq. IC.l 



^4 + 1) 

2 cosh 



J 4 = 



7RJ 4 (cj| + 4) 

6sinh(2ft) ■ 



(Cll) 



/i(l) = 
J X (3) = 



-2(5 

7 

-4J \ 
37 ) 



1/2 



C0Sh (^)' 

3/2 ^(9w 2 - 6) 



Suj cosh 



/ 3ttu \ 



h(2) = 



-35 
~2j 



2ttuj z 



-^sinh(^j) ' 

8tt(4u; 2 - 5) 
3(-5)sinh(^) 



(C12) 
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